
# generate table D2

# load packages
library(rio) # load data
library(tidyverse) # data manipulation
library(stargazer) # generate tables
library(brglm2) # estimate penalized logit
library(lmtest) # cluster standard errors
library(sandwich) # cluster standard errors

# set working directory
setwd("~/replication_files/")

# load data for analysis
data_with_dictionary <- import("data/full_data.csv") %>%
  mutate(across(c(imf_program,field_discovery,grad_school_econ_usa,iso3c,year), as.factor)) # treat factors as factors

# model 1, table D2
logit1 <- glm(article_iv_promotes_governance ~ 
                grad_school_econ_usa + fdi_performance_index_lag +
                imf_program + price_crudeoil_lag + price_crudeoil_difference + 
                resource_rents_lag + log_gdp_per_capita_lag + gdp_growth_lag + field_discovery_lag + 
                polyarchy + left_executive + protest + year + iso3c, data = data_with_dictionary, 
              family = "binomial"(link = "logit"), method = "brglmFit", slowit = 0.8, maxit = 500)

# cluster standard errors
logit1_r <- logit1 %>%
  coeftest(vcovHC(logit1, type = 'HC0', cluster =  ~ iso3c))

# model 2, table D2
logit2 <- glm(article_iv_mentions_resources ~ 
                grad_school_econ_usa + fdi_performance_index_lag +
                imf_program + price_crudeoil_lag + price_crudeoil_difference + 
                resource_rents_lag + log_gdp_per_capita_lag + gdp_growth_lag + field_discovery_lag + 
                polyarchy + left_executive + protest + year + iso3c, data = data_with_dictionary, 
              family = "binomial"(link = "logit"), method = "brglmFit", slowit = 0.8, maxit = 500)

# cluster standard errors
logit2_r <- logit2 %>%
  coeftest(vcovHC(logit2, type = 'HC0', cluster =  ~ iso3c))

# model 3, table D2
poisson1 <- glm(resource_mentions_absolute ~ 
                  grad_school_econ_usa + fdi_performance_index_lag +
                  imf_program + price_crudeoil_lag + price_crudeoil_difference + 
                  resource_rents_lag + log_gdp_per_capita_lag + gdp_growth_lag + field_discovery_lag + 
                  polyarchy + left_executive + protest + year + iso3c, data = data_with_dictionary,
                family = "poisson"(link = "log"), method = "brglmFit")

# cluster standard errors
poisson1_r <- poisson1 %>%
  coeftest(vcovHC(poisson1, type = 'HC0', cluster =  ~ iso3c))

# generate table D2
stargazer(logit1_r,logit2_r,poisson1_r, # type = "text", 
          omit = c("iso3c","year"), dep.var.labels = c("Consultation Promotes Natural Resource Governance",
                                                       "Consultation Mentions Natural Resources",
                                                       "Natural Resource Term Frequency"),
          covariate.labels = c("Technocratic Finance Minister = 1", "FDI Performance Index",
                               "IMF Program = 1", "Crude Oil Price, t--1", "Crude Oil Price, Delta",
                               "Resource Rents, t--1", "Log GDP Per Capita, t--1", "GDP Growth, t--1", "Field Discovery = 1, t--1",
                               "Polyarchy", "Left Executive", "Protest Count"))

# retrieve fit stats
stargazer(logit1,logit2,poisson1, # type = "text", 
          omit = c("iso3c","year"), dep.var.labels = c("Consultation Promotes Natural Resource Governance",
                                                       "Consultation Mentions Natural Resources",
                                                       "Natural Resource Term Frequency"),
          covariate.labels = c("Technocratic Finance Minister = 1", "FDI Performance Index",
                               "IMF Program = 1", "Crude Oil Price, t--1", "Crude Oil Price, Delta",
                               "Resource Rents, t--1", "Log GDP Per Capita, t--1", "GDP Growth, t--1", "Field Discovery = 1, t--1",
                               "Polyarchy", "Left Executive", "Protest Count"))




